program main
USE SHARE
USE RNUND_INT
USE UMIAH_INT
! USE UVMID_INT
IMPLICIT NONE
INTEGER :: dumfnm,ik,ip,lp,i,j,errcod, STATUS, iostat,T,K,P,rednfm,NBS,ib,js,ic
INTEGER :: bktr(900),IPARAM(7),dumper, jm,dumnfm,NP,jp,dnsm,ljm,ja,NYP,IS,IE,it
INTEGER :: ir,finind
DOUBLE PRECISION, DIMENSION(6000,900,10) :: ochars,bchars
DOUBLE PRECISION, DIMENSION(6000,900) :: oweig,orit,bweig,brit,opw
DOUBLE PRECISION, DIMENSION(900,10) :: chrsum,sschr
DOUBLE PRECISION, DIMENSION(900) :: sumsz,opret,mpret,vpret
DOUBLE PRECISION :: rdum,momdum,btmdum,medum,sew,sow,svw,ssow,ssew,ssvw,sw,dumome,lmdum
DOUBLE PRECISION:: sumpar(10),sspar(10),sumc(10),ssumc(10)
DOUBLE PRECISION:: U,betdum,sbdum,FUN,FOU,sbsw,FX,rsddum,lszdum,BMdum,LYSDUM,afydum
DOUBLE PRECISION::para(900,10),EU(900),ttu, optu, lastu
DOUBLE PRECISION:: H(10,10),G(10),THETA(10)
DOUBLE PRECISION:: dbsumw,snwtp,avgnw
INTEGER, DIMENSION(6000) :: IRVEC
EXTERNAL FCN,GRAD,HESS
OPEN(unit= 14, file='/Users/chrislamoureux/BSV/22/data/BSVUSE22.dat', STATUS ='OLD',iostat= errcod)
OPEN(unit= 17, file='/Users/chrislamoureux/BSV/22/data/output/bootstrap/M62/g8/RWTS.dat', STATUS='replace', iostat= errcod)
OPEN(unit= 21, file='/Users/chrislamoureux/BSV/22/data/output/bootstrap/M62/g8/Rtheta.dat', STATUS='replace', iostat= errcod)
OPEN(unit= 22, file='/Users/chrislamoureux/BSV/22/data/output/bootstrap/M62/g8/Rrets.dat', STATUS='replace', iostat = errcod)
OPEN(unit= 23, file='/Users/chrislamoureux/BSV/22/data/output/bootstrap/M62/g8/Rnstks.csv', STATUS='replace', iostat = errcod)
OPEN(unit= 33, file='/Users/chrislamoureux/BSV/22/data/output/bootstrap/M62/g8/RISrets.csv', STATUS='replace', iostat = errcod)
CALL ERSET(3,0,-1)
zero=0.0d0
one=1.0d0
two=2.0d0
NP=1         ! number of parameters (i.e., length of theta vector)
gamma=8.0d0
T=744   ! number of months in total.  This does not include the first data collecting period.
!         660 is through 12/14 Draft on my webpage (8.15.22) uses T=708.
! T=636     ! Matches the original LZ draft  (ends in 12/12)
P=180 ! The next month (P+1) is the starting month for out-of-sample portfolios.
NYP=(T-P)/12
!
DO i=1,T
   bktr(i)=0
   DO j=1,NP
      chrsum(i,j)=zero
      sschr(i,j)=zero
   END DO
   sumsz(i)=zero
END DO
!
DO
   READ (14,2166,iostat=status) jm,rdum,momdum,BMdum,lszdum,betdum,lysdum,afydum,rsddum,finind,dumome
            2166 format (i3,1x,8(f12.6,1x),i1,1x,f17.4)
   IF(status /= 0)  EXIT
   IF (jm /= ljm) THEN
      bktr(jm)=1
   ELSE
      bktr(jm)=bktr(jm)+1
   END IF
   orit(bktr(jm),jm) = rdum    !  read one stock sequentially
   oweig(bktr(jm),jm)= dumome
   sumsz(jm)=sumsz(jm)+dumome
!    ochars(bktr(jm),jm,1) = momdum   ! 1
!    ochars(bktr(jm),jm,1) = BMdum    ! 2
!    ochars(bktr(jm),jm,1)  = lszdum  ! 3
!    ochars(bktr(jm),jm,1)  = betdum  ! 4
!**** Do NOT USE:    ochars(ktr(jm),jm,1)  = lysdum    ******
   ochars(bktr(jm),jm,1)  = afydum  ! 5
!    ochars(bktr(jm),jm,3)  = rsddum  ! 6
   DO j=1,NP
      chrsum(jm,j)=chrsum(jm,j)+ochars(bktr(jm),jm,j)
   END DO
   ljm=jm
END DO
DO i=1,T
   PRINT 1126, i,bktr(i); 1126 format (2(2x,i4))
   WRITE (23,1126) i,bktr(i)
END DO
! Standardize and normalize the parameters and define the market weights.
DO i=1,T
   DO j=1,NP
      chrsum(i,j)=chrsum(i,j)/bktr(i)
   END DO
   do j=1,bktr(i)
      oweig(j,i)=oweig(j,i)/sumsz(i)
      WRITE (17,1766) i,j,oweig(j,i); 1766 format (i3,2x,i4,f12.9)
   end do
END DO
DO i=1,T
   DO j=1,bktr(i)
      DO ja=1,NP
         sschr(i,ja)=sschr(i,ja)+(ochars(j,i,ja)-chrsum(i,ja))*(ochars(j,i,ja)-chrsum(i,ja))
      END DO
   END DO
END DO
DO i=1,T
   DO j=1,NP
      sschr(i,j)=sqrt(sschr(i,j)/(bktr(i)-one))
   END DO
END DO
DO i=1,T
   DO j=1,bktr(i)
      DO ja=1,NP
         ochars(j,i,ja)=(ochars(j,i,ja)-chrsum(i,ja))/sschr(i,ja)
      END DO
   END DO
END DO
!
!   out of sample series
!   IF (T /= 100) STOP 1091
NBS=10000
DO ib=1,NBS
! Construct the pseudo sample for this bootstrap:
   DO jm=1,T
      DO ic=1,NP
         sumc(ic)=zero
         ssumc(ic)=zero
      END DO
      sbsw=zero
      CALL RNUND(bktr(jm),irvec)
      DO js=1,bktr(jm)
         brit(js,jm)=orit(irvec(js),jm)
         bweig(js,jm)=oweig(irvec(js),jm)
         sbsw=sbsw+bweig(js,jm)
         DO ic=1,NP
            bchars(js,jm,ic)=ochars(irvec(js),jm,ic)
            sumc(ic)=sumc(ic)+bchars(js,jm,ic)
         END DO
      END DO
      DO ic=1,NP
         sumc(ic)=sumc(ic)/bktr(jm)
      END DO
      DO js=1,bktr(jm)
         bweig(js,jm)=bweig(js,jm)/sbsw
         DO ic=1,NP
            ssumc(ic)=ssumc(ic)+((bchars(js,jm,ic)-sumc(ic))*(bchars(js,jm,ic)-sumc(ic)))
         END DO
      END DO
      DO ic=1,NP
         ssumc(ic)=sqrt(ssumc(ic)/bktr(jm))
      END DO
      DO js=1,bktr(jm)
         DO ic=1,NP
            bchars(js,jm,ic)=(bchars(js,jm,ic)-sumc(ic))/ssumc(ic)
         END DO
      END DO
   END DO
!
! At this point a bootstrap sample is set up.
!
   DO jp=1,NP
     theta(jp)=zero
     sumpar(jp)=zero
     sspar(jp)=zero
   END DO
   DO k=1,NYP           ! We update theta at the beginning of each year using historical data starting at 
      IS=(k-1)*12+1     !    the beginning of the rolling sample.
      NMOS=P            ! NMOS is the length of the in-sample period (constant when we roll).
      IE=P+(k-1)*12     ! End of the rolling period.
      it=0
      DO ir=IS,IE
         it=it+1
         ktr(it)=bktr(ir)                          ! Even though the bootstrap sample is set up I have to construct
         DO js=1,bktr(ir)                          ! an additional data set whose month 1 corresponds to lagging 180
            rit(js,it)=brit(js,ir)                 ! months from the start of the oos period.
            weig(js,it)=bweig(js,ir)               !  ==> I have to change the ir subscript to it for the following
            DO ic=1,NP                             !      variables: ktr, iand for each stock: rit, weig, and chars
               chars(js,it,ic)=bchars(js,ir,ic)
            END DO
         END DO
      END DO
!
!       PRINT 4266, NMOS,P; 4266 format (' in main before umiah nmos: ',i4,' p ',i4)
!
      IPARAM(1)=0
      CALL UMIAH (FCN,GRAD,HESS,theta,N=NP,IPARAM=IPARAM,FVALUE=FOU)   ! IMSL Math V2 p. 1213.
      WRITE (33,3366) k; 3366 format (i4)
      DO i=1,NMOS
         WRITE (33,3367) i,ISOPRT(i); 3367 format (i4,1x,f14.10)
      END DO
!       PRINT 1426, FOU,(theta(jp),jp=1,NP),IPARAM(3)
            1426 format (' UMIAH finished Max U: ',f12.7,' theta: ',3f12.7,' no iterations: ',i4)
           FOU=FX
           EU(k)=FOU
           WRITE(21,114) k, (theta(jp),jp=1,NP); 114 FORMAT(1X, I8, 10(2X, F12.8))
!             WRITE(21,114) k, (theta(jp),jp=1,NP); 114 FORMAT(1X, I8, 10(2X, F14.6))
         DO i = 1, 12   ! assign parameters to the next year
            DO jp=1,NP
               para(IE+i,jp) = theta(jp)
               sumpar(jp)=sumpar(jp)+ theta(jp)
            END DO
         END DO
   END DO   ! Loop over the out-of-sample years
!
!    Optmizations complete.  Set up the returns for this bootstrap sample.
!
   U=zero
   avgnw=zero
   DO i=P+1,T
      snwtp=zero
      sw=zero
      DO j = 1, bktr(i) ! optimal weights and returns of stocks
         sc=zero
         DO jp=1,NP
            sc = sc+para(i,jp)*bchars(j,i,jp)
         END DO
         opw(j,i)=bweig(j,i)+sc/bktr(i)
         IF (opw(j,i) < zero) THEN
            snwtp=snwtp+opw(j,i)
         END IF
         sw=sw+opw(j,i)
!        WRITE(20,121) MDATE(i),bktr(i),weig(j,i),opw(j,i); 121 FORMAT(1X,I8,1X, I5,2(2X,F16.13))
      END DO
      avgnw=avgnw+snwtp
!  PRINT 5632, i,MDATE(i),sw; 5632 format (' Month ',i4,' date: ',i8,' sum wts: ',f12.10)
      opret(i) = zero; mpret(i)=zero; vpret(i)=zero
      dbsumw=zero
      DO j =1,bktr(i) ! Portfolio returns in each month
         opret(i) = opret(i) + opw(j,i)*brit(j,i)
         mpret(i) = mpret(i) + brit(j,i)
         vpret(i) = vpret(i) + bweig(j,i)*brit(j,i)
         dbsumw=dbsumw+opw(j,i)
      END DO
!       PRINT 8761, i,dbsumw; 8761 format (2x,i4,' debugging sum of opt wts: ',f12.9)
      U = U +((one+opret(i))**(one-gamma))/(one-gamma)
   END DO
   U=U/(T-P)
   avgnw=avgnw/(T-P)
   PRINT 5439, U; 5439 format (' Mean oos utility: ',f12.9)
   PRINT 5446, avgnw; 5446 format (' Avg neg wts: ',g16.7)
   WRITE (17,1776) avgnw; 1776 format (f16.8)
   sew=zero
   svw=zero
   sow=zero
   ssew=zero
   ssvw=zero
   ssow=zero
!    DO i=1,T
!       PRINT 7651, i,MDATE(i),bktr(i); 7651 format (1x,i4,2x,i8,2x,i7)
!    END DO
   DO i = P+1,T
      mpret(i)=mpret(i)/bktr(i)
      sew=sew+mpret(i)
      svw=svw+vpret(i)
      sow=sow+opret(i)
      WRITE(22,117) i,opret(i),mpret(i),vpret(i); 117 FORMAT(1X,I8,4(2X,F12.8))
   END DO
   sew=sew/(T-P)
   svw=svw/(T-P)
   sow=sow/(T-P)
   DO i=1,NP
      SUMPAR(i)=SUMPAR(i)/(T-P)
   END DO
   DO i=P+1,T
      ssew=ssew+(mpret(i)-sew)*(mpret(i)-sew)
      ssvw=ssvw+(vpret(i)-svw)*(vpret(i)-svw)
      ssow=ssow+(opret(i)-sow)*(opret(i)-sow)
      DO jp=1,NP
         sspar(jp)=sspar(jp)+(para(i,jp)-sumpar(jp))*(para(i,jp)-sumpar(jp))
      END DO
   END DO
   DO jp=1,NP
      sspar(jp)=sqrt(sspar(jp)/(T-P))
   END DO
   ssew=sqrt(ssew/(T-P))
   ssvw=sqrt(ssvw/(T-P))
   ssow=sqrt(ssow/(T-P))
!   PRINT 7861, gamma; 7861 format (' Gamma ',f8.3)
   PRINT 6574, ib; 6574 format (' Bootstrap sample: ',i6)
   PRINT 7865; 7865 format (28x,' Mean ',8x,'Standard Deviation')
   PRINT 7866, sew,ssew; 7866 format(' Equally weighted Index',2(4x,f12.9) )
   PRINT 7867, svw,ssvw; 7867 format(' Value weighted Index  ',2(4x,f12.9) )
   PRINT 7868, sow,ssow; 7868 format(' BSV Optimal Portfolio ',2(4x,f12.9) )
   DO jp=1,NP
      PRINT 9436, jp,sumpar(jp),sspar(jp); 9436 format (1x,i2,2(2x,f13.8))
   END DO
! 1 bootstrap sample finished.
END DO
STOP; END
